rm(list = ls())
cat("\014")
library(tidyverse, quietly = T)
library(sf, quietly = T)
library(gganimate, quietly = T)
library(openair, quietly = T)
library(lubridate, quietly = T)
library(zoo, quietly = T)
library(ggthemes, quietly = T)
library(ggmap, quietly = T)
library(sp, quietly = T)
library(magick, quietly = T)
Import the time-activity diary for the line and tidy it up.
diary <- read_csv('air_quality_data/victoria_line_2014-12-08.csv', col_types = cols())
base <- tibble(date_time = seq(min(diary$start_time), max(diary$end_time), by = 'min'))
base <- left_join(base, diary, by = c('date_time' = 'start_time')) %>%
select(-fake_id)
base <- left_join(select(base, -end_time),
filter(base, date_time != end_time),
by = c('date_time' = 'end_time')) %>%
mutate(line.x = if_else(is.na(line.x) & !is.na(line.y), line.y, line.x),
station.x = if_else(is.na(station.x) & !is.na(station.y), station.y, station.x),
environment.x = if_else(is.na(environment.x) & !is.na(environment.y), environment.y, environment.x),
sub_environment.x = if_else(is.na(sub_environment.x) & !is.na(sub_environment.y), sub_environment.y, sub_environment.x)) %>%
select(date_time, line.x, station.x, environment.x, sub_environment.x) %>%
rename(line = line.x, station = station.x, environment = environment.x, sub_environment = sub_environment.x)
rm(diary)
Import the measured air quality to join to the time-activity diary
air_quality <- read_csv('air_quality_data/LUexportMay2016.csv', col_types = cols()) %>%
rename_all(tolower) %>%
mutate(datetime = as.POSIXct(datetime, format='%d/%m/%Y %H:%M')) %>%
filter(sitecode == 'CAR')
Join the time activity and air quality
line <- left_join(base, air_quality, by = c('date_time' = 'datetime')) %>%
mutate(scaledvalue = if_else(species == 'PM25', value, scaledvalue)) %>%
select(-value) %>%
rename(concentration = scaledvalue)
rm(air_quality, base)
Import the background air quality for that day to correct the PM2.5 measurements
background_pm25 <- importKCL(site = "kc1", year = c(2014,2015,2016), pollutant = "pm25", met = FALSE,
units = "mass", extra = FALSE)
## NOTE - mass units are used
## ug/m3 for NOx, NO2, SO2, O3; mg/m3 for CO
## PM10_raw is raw data multiplied by 1.3
background_pm25 <- data.frame(background_pm25, day = as.Date(format(background_pm25$date)))
background_pm25 <- aggregate(pm25 ~ day, background_pm25, mean)
Join the background PM2.5 to the monitored concentrations
line <- mutate(line, day = date(date_time)) %>%
left_join(background_pm25, by = c('day' = 'day')) %>%
select(-day)
Correct the PM2.5 data
line <- mutate(line, concentration =
case_when(concentration > pm25/0.44 & species == 'PM25' & !is.na(concentration) ~ pm25 + (1.82 * (concentration - (pm25/0.44))),
concentration <= pm25/0.44 & species == 'PM25' & !is.na(concentration) ~ concentration * 0.44,
TRUE ~ concentration)) %>%
select(-pm25)
Get the station locations
stations <- st_read('https://raw.githubusercontent.com/dracos/underground-live-map/master/bin/stations.kml', quiet = T) %>%
select(-Description) %>%
rename_all(tolower) %>%
mutate(name = gsub(' Station', '', name)) %>%
mutate(name = gsub("'", '', name)) %>%
st_zm(drop = TRUE)
Join locations to the concentrations and station names
line <-left_join(line, stations, by = c('station' = 'name')) %>% st_as_sf()
Fill in the missing coordinates using linear interpolation
st_geometry(line) <- as_tibble(st_coordinates(line)) %>%
mutate(X = na.approx(X), Y = na.approx(Y)) %>%
st_as_sf(coords = c('X', 'Y')) %>%
st_geometry()
Start trying to animate. First the graph animation.
text_placement <- filter(line,species == 'PM25') %>%
as_tibble() %>%
select(-geometry) %>%
group_by(species) %>%
summarise(max_x = max(date_time), max_y = max(concentration))
line <- left_join(line, text_placement, by = c('species' = 'species'))
pm25 <- filter(line,species == 'PM25' & !is.na(concentration))
pm25$moving_mean <- NA
pm25$pc25 <- NA
pm25$pc75 <- NA
for (i in 2:nrow(pm25)) {
pm25[i,'moving_mean'] <- as.integer(round(mean(pm25[1:i,]$concentration),0))
pm25[i,'pc25'] <- as.integer(round(quantile(pm25[1:i,]$concentration, 0.25),0))
pm25[i,'pc75'] <- as.integer(round(quantile(pm25[1:i,]$concentration, 0.75),0))
}
graph <- ggplot(data = filter(line,species == 'PM25')) +
geom_path(aes(date_time, concentration), colour = '#0099CC', size = 1.2) +
geom_ribbon(data = filter(pm25,!is.na(pc25)), aes(x=date_time, ymin = pc25, ymax = pc75), alpha = 0.2) +
geom_hline(data = filter(pm25,!is.na(pc25)), aes(yintercept = moving_mean), alpha = 0.2, colour = 'red') +
geom_point(aes(date_time, concentration), colour = 'red', size = 4) +
geom_text(data = filter(pm25,!is.na(pc25)), aes(x = max(date_time), y = moving_mean, label = moving_mean), hjust = 0, size = 8, colour = 'red') +
geom_vline(xintercept = as.POSIXct('2014-12-08 10:00'), size = 1) +
coord_cartesian(clip = 'off') +
geom_text(aes(min(date_time), max_y,label = round(concentration,0)),
size = 10, hjust = 0, vjust = 1, colour = '#0099CC') +
theme(axis.text.x = element_blank(),
axis.text.y = element_text(size = 12),
axis.title.y = element_text(size = 12),
axis.line = element_line(colour = 'black'),
plot.margin = margin(5.5, 40, 5.5, 5.5)) +
xlab('') +
ylab('PM2.5 ug/m3') +
transition_reveal(date_time)
graph_animation <- animate(graph, fps = 3, height = 600, renderer = gifski_renderer(loop = FALSE))
anim_save(file='graph_animation.gif', animation = graph_animation, path = 'outputs')
graph_animation
Now the spatial animation.
line <- st_set_crs(line,4326)
map_area <- c(as.vector(st_bbox(line))[1]-0.01, as.vector(st_bbox(line))[2]-0.01, as.vector(st_bbox(line))[3]+0.01, as.vector(st_bbox(line))[4]+0.01)
map <- ggmap(get_map(location = map_area), source = "google", maptype = "roadmap", zoom = 13) +
geom_path(data = as_tibble(st_coordinates(filter(line,species == 'PM25')[1:56,])), aes(X, Y), colour = "#0099CC", size = 1) +
geom_sf(data = filter(line,species == 'PM25'), size=4, colour = 'red',inherit.aes = FALSE) +
coord_sf(datum=NA) +
transition_time(date_time) +
ease_aes('linear') +
theme(axis.title = element_blank(),
panel.background = element_blank(),
panel.border = element_rect(colour = 'black', fill = NA),
axis.text = element_blank())
map_animation <- animate(map, fps = 3, height = (470*2), width = (200*2))
anim_save(file='map_animation.gif', animation = map_animation, path = 'outputs')
map_animation